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' Abstract. I discuss the conformal approach to the numerical simulation of radiating 

, isolated systems in general relativity. The method is based on conformal compactifica- 

tion and a reformulation of the Einstein equations in terms of rescaled variables, the 
?H ' so-called "conformal field equations" developed by Friedrich. These equations allow to 

include "infinity" on a finite grid, solving regular equations, whose solutions give rise 
to solutions of the Einstein equations of (vacuum) general relativity. The conformal 
approach promises certain advantages, in particular with respect to the treatment of 
radiation extraction and boundary conditions. I will discuss the essential features of 
the analytical approach to the problem, previous work on the problem - in particu- 
lar a code for simulations in 3+1 dimensions, some new results, open problems and 

■ strategies for future work. 

t> 
in 

1 Introduction 

' In order to understand the physical content of the theory of general relativity, it 

is desirable to both mathematically understand its solutions and observationally 
""^ ' understand the physical phenomena for which the theory is relevant. The latter 

effort typically requires predictions from the theory, both qualitative and quan- 
^ ' titative - such as gravitational wave templates or binary pulsar deceleration 

(3JT), parameters. The lack of genericity in available exact solutions then naturally 

leads to the use of approximation methods such as post-Newtonian approxi- 
mations, perturbation theory or numerical analysis, which allows very general 

■ non-perturbative approximations. Concrete solutions do however also play an 
^ , important role in the quest for a mathematical understanding of the solution 

space. The experience gained from such solutions can suggest theorems, test 
conjectures, or lead to the discovery of previously unknown phenomena. For 
some particularly interesting examples see 0] or ||. The construction and 
study of solutions, be it with approximate or exact methods, obviously profits 
from a sound mathematical basis in the form of well-posed equations, analytic 
estimates and the likes. Eventually - hopefully - it will also profit from obser- 
vational evidence! 

In the following I will discuss a particular approach to the numerical solu- 
tion of the Einstein field equations, which addresses the problems associated 
with the treatment of asymptotic regions by conformal compactification. The 
interest in asymptotic regions is rooted in the problem of describing isolated 
systems. Physical intuition suggests that many astrophysical processes (whether 
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they are of actual astrophysical relevance or rather hypothetical) should essen- 
tially be independent of the large-scale structure of the universe, or, say, the 
local galaxy. The idealization of an isolated system, where the geometry ap- 
proaches a Minkowski geometry at large distances, thus forms the basis for the 
general-relativistic analysis of processes which are essentially of non-cosmological 
nature. The mathematical formalization of the physical idea of isolated systems 
is the concept of asymptotically flat spacetimes. This formalization is already 
nontrivial, due to the lack of a preferred background geometry or coordinate 
system - with respect to which one could define "distance" and the appropri- 
ate limits. Conformal compactification, however, renders possible a discussion 
of asymptotically flat spacetimes in terms of local differential geometry. In this 
approach, pioneered by Penrose ||, an unphysical Lorentzian metric g ab is in- 
troduced on an unphysical manifold M. which gives rise to the physical metric 
gab by the rescaling g ab = f2~ 2 g ab . The physical manifold M is then given by 
M = {p G M. | f2(p) > 0}. In this picture physical "infinity" corresponds to 
a three-dimensional boundary of a four-dimensional region in A4, defined by 
Q = Q. Limiting procedures and approximations can thus be replaced by local 
differential geometry on the boundary. 

In gravitational theory, quantities such as the total mass, (angular) momen- 
tum or emitted gravitational radiation can only consistently be defined at "in- 
finity" . In the conformal approach the unambiguous extraction of gravitational 
waves from a numerical spacetime is straightforward. In the "traditional" ap- 
proach to dealing with asymptotic falloff in numerical relativity, where one in- 
troduces an arbitrary spatial cutoff, matters are much more complicated and 
ambiguities are introduced which one would have to get rid off by complicated 
limiting procedures. Without at least being able to define a clean concept of radi- 
ation leaving or entering a system, it is furthermore very hard to define physically 
realistic and consistent boundary conditions at finite distance. The traditional 
approach is thus not completely satisfactory both from a mathematical but also 
from a practical point of view. Here we discuss the principal ideas of the envi- 
sioned "conformal cure" , the technical and conceptual problems associated with 
it, and the current status of this approach. 

It is easy to see that the conformal cure can not be straightforward, by writing 
Einstein's vacuum equations in terms of Q and g ab '- 

Gab[^ 2 g a b] = G ab [g ab ] + 1 (V a V 6 X2 + <? afc V c V c l2) + g ab (V /2) V c f2 . (1) 

This expression is singular for Q = 0, multiplication by Q 2 also does not help 
here because then the principal part of the partial differential equations encoded 
in G a b would degenerate at fl = 0. The conformal compactification approach 
thus can not be carried to the level of the field equations in a straightforward way. 
This step however has been achieved by Friedrich, who has developed a judicious 
reformulation of the equations |^J^ , p?|p7[|2c| . These conformal field equations 
are regular equations for g ab and certain additional independent variables. 

In analytical work, such global methods have proven to provide essential 
simplifications leading to new results and insights. Already by providing a dif- 
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ferent point of view on some of the essential problems in numerical relativity, 
the conformal picture is quite helpful and can stimulate new ideas. Certainly, we 
desire more - to make this approach also a practical tool. There is significant 
hope, that global methods will eventually show advantages for practical numeri- 
cal work, and despite the small number of researchers involved so far (may there 
be more!), some significant progress in this direction has been made. 

In the present article I will try to sketch the present status of the quest for the 
conformal cure and discuss some important open questions. We will start with 
a brief introduction of the concepts of asymptotic flatness in terms of conformal 
compactification in Sec. highlighting some important features of "future null 
infinity" , and then discuss the conformal field equations. In Sec. || I will discuss 
some explicit examples of compactifying Minkowski spacetime, both to paint a 
more concrete picture of our scenario, and to set the arena for some numerical 
code tests. Sec. |] contains a brief overview of the history of numerical work on 
the conformal field equations, leading to a description of a 3D code written by 
Hubner f||6|Q|| . New results from 3D calculations performed with this code will 
be presented in Sec. [|, and a discussion will be given in Sec. |^, concluding with 
a roadmap for future work. 



2 Compactification and the Mathematical Description of 
Isolated Systems 



The material in this section is intended to present some essential ideas in a 
condensed form. The reader should be aware that I am not doing justice here to 
subtleties and long history of the mathematical description of isolated systems 
in general relativity - rather this section intends to motivate to look into more 
complete reviews such as PHQ. 



2.1 Asymptotic Flatness and Compactification 

As noted above the formulation of the concept of asymptotic flatness is far from 
straightforward in GR, due to the absence of a background metric or preferred 
coordinate system, in terms of which falloff rates can be specified. A resolution 
of this problem is provided by a definition of asymptotic flatness, where, after 
a suitable conformal rescaling of the metric, "points at infinity" are added to 
the manifold. One thus works on a compactified auxiliary manifold, and local 
differential geometry can be used to study the asymptotic properties of the 
gravitational field. We will give a simple definition of asymptotic flatness here, 
which for our purposes catches all essential features. For alternative definitions 
and more detailed explanations compare for example [§0 ^0 , . 

Definition 1 (asymptotic simplicity) 

A smooth spacetime (M , g a b) is called asymptotically simple, if there exist an- 
other smooth spacetime (A4,g a b) and a scalar function [2 such that: 
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1. M is an open submanifold of M. with smooth boundary 8A4 — J? (Scri). 

2. g ab = [2 2 g a b on M, with Q > on M, Q = on J and V ' a fl ^ on J . 

3. Every null geodesic in M. acquires two end points on J? . 

Definition 2 (asymptotic flatness) 

Asymptotically simple spacetimes are called asymptotically flat if their Ricci 
tensor R a b vanishes in a neighborhood of J? . 

Examples of asymptotically simple spacetimes, which are not asymptotically 
flat are the de Sitter and anti-de Sitter solutions. Correspondingly to asymptoti- 
cally flat spacetimes one can consider asymptotically de Sitter and anti-de Sitter 
spacetimes. Note that the completeness condition 3 in Def. [l], which ensures that 
the entire boundary is included, excludes black-hole spacetimes. For modifica- 
tions to weaken condition 3, thus allowing black holes, see the definitions of E3| 
or |1. For example, the definition of weak asymptotic simplicity fl^j requires 
condition 3 to hold only in a neighborhood of J? . See e.g. 0j for a discussion of 
asymptotic flatness at spacelike infinity (i.e. the part of infinity which is reached 
along spacelike geodesies) versus null infinity (i.e. the part of infinity which is 
reached along null curves). The notion of asymptotic flatness at timelike infinity 
does not make much sense in a general situation, because then all energy would 
have to be radiated away, leaving only flat space behind - excluding black holes 
or "stars" . For weak data however, in vacuum say, where all radiation eventually 
disperses, once expects asymptotic flatness to hold also at timelike infinity, this 
issue will be discussed below in application to concrete spacetimes. 

The notion of asymptotic flatness of isolated systems turns out to be inti- 
mately related to the possibility of defining the total energy-momentum for such 
systems in general relativity - remember that no well-defined local energy den- 
sity of the gravitational field is known (compare e.g. Sec. 11.2 of the textbook 
of Wald |^]). However, total energy-momentum quantities, which transform as a 
4- vector under asymptotic Lorentz transformations, can be assigned to null and 
spatial infinity of asymptotically flat spacetimes. If a manifold has more than 
one asymptotically flat end, e.g. in the presence of wormholes of the Einstein- 
Rosen-bridge type, then different energy-momenta can be associated with each 
of these asymptotic regions. 

The expression for the energy momentum four-vector at spatial infinity has 
been given first by Arnowitt, Deser and Misner in 1962 Jll| in the context of 
the Hamiltonian formalism, and is usually called the ADM momentum, the time 
component being called ADM mass. The ADM energy corresponds to the energy 
of some Cauchy surface, i.e. a snapshot of the spacetime at some fixed time. It is 
a constant of motion and can therefore be expressed in terms of the initial data 
on an asymptotically flat Cauchy hypersurface. 

The expression for the energy momentum at null infinity, usually referred to 
as the Bondi energy-momentum, can be associated with a fixed retarded time, 
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i.e. some asymptotically null surface. The decrease of this quantity measures 
the energy-momentum carried away by gravitational radiation. For a brief in- 
troduction and references to original work on different definitions of the Bondi 
mass see e.g. the textbook of Wald The formulation most appropriate for 
usage in numerical codes based on the conformal field equations was given by 
Penrose j4| , and defines the Bondi mass in terms of the behavior of certain pro- 
jections of the Weyl tensor at and the shear of the outgoing congruence of 
null geodesies orthogonal to in the gauge defined below by (||) . It was already 
shown in 1962 by Bondi, van der Burg and Metzner JlJ] that the Bondi mass 
Mb can only decrease with time: gravitational radiation always carries positive 
energy away from a radiating system. Note that this means in particular, that 
while compactification at spatial infinity would lead to a "piling up" of waves, 
at J? + this effect does not appear. In the compactified picture the waves leave 
the physical spacetime through the boundary J> + . 

A fundamental issue of general relativity is the positivity of the ADM and 
Bondi energies. Although it is trivial to write down a metric with negative mass 
if no conditions on the energy momentum tensor are imposed, for reasonable 
matter fields with nonnegative energy density (thus satisfying the dominant 
energy condition), non-negativity of the ADM and Bondi energies is expected 
on physical grounds: if the energy of an isolated system could be negative, it 
would most likely be unstable and decay to lower and lower energies. Indeed, a 
proof of the positive definiteness of the ADM energy has been given in 1979 by 
Schoen and Yau |ll| (several simplified proofs have been given later), and was 
extended to the Bondi mass in 1982 by Horowitz and Perry f[6| |. 

2.2 What is Scri? 

We will now have a closer look at J 1 and discuss some of its features, which 
will allow us to understand the basic ideas of radiation extraction and help us to 
understand some issues related with choosing boundary conditions for numerical 
solutions of the conformal field equations. 

Looking at (|l|) and multiplying by Q 2 , one can see that for a vacuum space- 
time, G a b — 0, (V c i7)V c i7 = at J?, which thus must consist of null surfaces. 
In fact, one can then prove (see e.g. 0), that 

1. J 1 has two connected components, each with topology S 2 x R. 

2. The connected components of are smooth null hypersurfaces in A4, and 
as such are generated by null geodesies. 

3. The congruence of null geodesic generators of is shear free. 

The two connected components are called future null infinity and past 

null infinity («^ _ ), and provide the future and past endpoints for null geodesies 
in M.. In a naive picture they could be viewed as emanating from a point i° 
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which represents spatial infinity [j. These features will become more graphic 
when dealing with explicit examples below. 

Note that there is gauge freedom in the choice of the conformal factor: one 
is free to rescale the conformal factor fl by some u> > such that Q = tuft, 
g a b — LU 2 g a b — f2 2 g a b- It is an interesting exercise (see Sec. 11.1 of to prove 
that outside any neighborhood of i° - on say - one can always use this 
conformal gauge freedom to achieve 

V a V 6 r) = on/+, (2) 

where V a is the derivative operator compatible with the metric g a b- This con- 
formal gauge implies, that the null tangent n a = g ab Vb& to the null geodesic 
generators of satisfies the affincly parameterized geodesic equation, 

n a V a h b = . (3) 

Consequently, expansion of the generators of ^ vanishes in addition to the shear 
and twist (n a is a gradient). Using the remaining gauge freedom of u>, we can 
choose coordinates such that the metric on Jt takes the form 

ds 2 U+ = + d9 2 + sin 2 ddcf) 2 , (4) 

where u is the affine parameter of the null geodesic generators, scaled such that 
n a V a M = 1 (see e.g. Chpt. 11 of ||). The cuts f] of of constant u thus become 
metric spheres. The coordinate u is generally known as Bondi parameter or 
Bondi time. The conformal gauge ^ and the coordinates (^) prove very useful 
in the analysis of the geometry in a neighborhood of J* - in particular for the 
extraction of radiation. The existence of a natural time coordinate (at least up 
to affine transformations along each generator) is very interesting for numerical 
applications, where at least asymptotically one can get rid of much of the slicing 
arbitrariness of the interior region. It is nontrivial but rather straightforward 
to actually (numerically) find this gauge of J^ + , which is also required by the 
standard formulas to compute the energy-momentum at ^ + and the emitted 
radiation - to be given below. 

Before discussing how to compute the radiation, it is useful to idealize a de- 
tector (here I will follow the discussion in [|lj ) . In physical space - far away form 
the sources - we could think of a detector as a triad of spacelike unit vectors 
attached to the worldline of some (timelike) observer. Let us further assume for 
simplicity that the observer moves along a timelike geodesic parametrized by 
proper time and that the triad is transported by Fermi- Walker transport. It is 
not hard to show - see Frauendiener pi] , that taking the appropriate limit in 
the compactified spacetime, the observer worldline converges to a null geodesic 

1 The structure of i is however quite subtle, significant progress toward its under- 
standing in terms of the field equations has recently been achieved by Friedrich 

2 A cut of J 2 " is a two-dimensional spacelike cross section of which meets every 
generator once. 
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generator of J^ + . Taking the limit along a Cauchy surface it converges to the 
point i , where one could naively expect an observer to end up when shifted to 
larger and larger distances (this limit is however not appropriate in the context 
of computing the radiation). Furthermore, the proper time parameter of the 
observer converges to Bondi time. The arbitrariness of boosting the observers 
is reflected in the affine freedom of choosing the Bondi parameter at <P . The 
description of ^ + thus could be condensed into the statement that it idealizes 
us - the observers of astrophysical phenomena happening far away. By working 
with the idealization, the approximations and ambiguities associated with de- 
tectors at a finite distance have transformed into a surprisingly simple geometric 
picture! Note that this simplification has to be taken with the typical care re- 
quired in the treatment of idealizations in (theoretical) physics: Under practical 
circumstances, e.g. computing the actual signal at a gravitational wave detector, 
J 2 " more realistically corresponds to an observer that is sufficiently far way from 
the source to treat the radiation linearly, but not so far away that cosmological 
effects have to be taken into account. In order to compute the detected signal in a 
realistic application, cosmological data and the fact that an earthbound detector 
moves in a complicated way relative to the source all have to be considered. 

We will next discuss a "detector-frame" adapted to - the commonly used 
Bondi frame. For a much more complete discussion of Bondi-systems see e.g. the 
excellent review by Newman and Tod JToj ] . There a characteristic framework is 
used to set up the Bondi frame in a whole neighborhood of which is nec- 
essary to compute derivatives, entering e.g. the definition of the spin coefficient 
a defined below in (EJ). In the current approach, the Bondi system is only de- 
fined at J^ + : initial data can be set up, such that all necessary quantities can 
be propagated along the generators of (2^] . 

With being a null surface, it is most natural to use a null frame, consisting 
of 2 null vectors and 2 spacelike vectors x a , y a , which can be considered as the 
idealizations of the arms of an interferometric gravitational wave detector. The 
vectors x a , y a are commonly treated in the form of two complex null vectors 
m a , rh a , with 

a a i ■ a a — 1 

m = x + ty , m m a = 1 , 

where x a and y a are real vectors tangent to the cuts of The null vectors 
are taken as the affine tangent n a and l a = V a it, which satisfy 



The tetrad vectors l a , m a and fh a are parallely propagated along the generators, 
which yields transport equations that define them on all of J? + once initial values 
are chosen. 

The Bondi-mass can then be computed in terms of the spin-coefficient a and 
the rescaled Weyl tensor components -02 and ^4: 

a = g ab l a m c \7 c m b , (5) 
^2 = d abcd l a m b fh c n d , (6) 
V^4 = d abcd n a m b n c rh d . (7) 
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In terms of these quantities the Bondi mass can be defined as 

Mb = — 7=3 / (^2 + era) dA , (8) 

V47T J 

the outgoing radiation can be computed to be 

M B = — ^=3 / {&a) dA , (9) 

V47T J 

where A is the area of the cuts of J?^ ', and / = n a V a / = d u f. Furthermore, 

b = -$ 4 

can be used to evolve a, where both a and a can be computed on the initial 
slice. 

This procedure has been implemented by Hiibner and Weaver p3 for 2D 
codes and the 3D code used to obtain the results in Sec. ||, and has been tested 
and proven accurate for several types of spacetimes ^] . Frauendiener describes 
his implementation and some results in jl9). There are two essential problems 
in these implementations: First of all, the gauge conditions will not usually re- 
sult in a slicing of by cuts of constant Bondi time u. This means that 
interpolation has to be used between different slices of the numerical evolution. 
Second, in those formulations of the conformal field equations that have so far 
been used in numerical implementations, the conformal factor fl is an evolution 
variable and not specified a priori, J? will in general not be aligned with grid 
points. This results in further technical complications and an additional need 
for interpolation. When dealing with the physically interesting case of a of 
spherical topology, at least two patches have to be used to represent the Bondi 
tetrad (l a ,n a ,m a ,fh a ). Frauendiener has achieved to control the movement of 
J r+ through the grid by the gauge choice for his formulation Q, in particu- 
lar the shift vector can be chosen such that J 1 does not change its coordinate 
location. 



2.3 The Conformal Field Equations 

Several formulations of the conformal field equations are available, the main 
difference being whether the conformal factor fl can be specified a priori or is 
determined as a variable by the equations. In the original formulation |2^ , p5| 
and its descendants ||f]||[l8| fl (and derivatives) are evolved as dependent vari- 
ables. All existent numerical codes are based on equations of this type. A later 
version of the equations allows to fix fl a priori, and has been used to develop 
a new treatment of spatial infinity i° [ p7|j2^ , |29|| . However, the formulation and 
treatment of these equations is more involved, and its numerical solution has not 
yet been attempted. 

In the following we discuss a metric based formulation of the "original" ver- 
sion of the conformal field equations, which forms the basis for Hubner's codes 

mum- 
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When deriving the conformal field equations, it turns out to be useful to start 
with the splitting of the Riemann tensor into its trace-free (the Weyl tensor) and 
trace (Ricci tensor and scalar) parts. Additionally we define the tracefree Ricci 
tensor R a b = R ao — \ g a b R and the rescaled Weyl tensor 

d a bc d = ft 1 C a b c d ■ (10) 
The requirement that the physical scalar curvature R vanishes implies 

6nv a v a f2 = i2(v a n) (v a n) - n 2 R , (n) 

Note that this equation is not manifestly regular at Q = 0, but it is actually 
possible to show that if (pi]) is satisfied at one point, then by virtue of the other 



equations (12,15 , 1 ^ 1_7|, 1 8h t o be given below, it has to be satisfied everywhere. 
The whole system~( |l 1| Jl 2| , |1 | , |14j , |l7| , |l8| ) is then regular in the sense that this point 
does not have to be located at . The vacuum Einstein equations R a b = 
then yield 

V a Vbft = X -g ab V c V c ft - i Rab ft ■ (12) 
Finally, commuting covariant derivatives in the expression 

g bc V c V b V a ft 

and then using (n2h again yields 



^V a (V b Vbft) =~R ab V b ft~^ftVaR- -^VaftR- (13) 
Equations for the metric can be obtained by the identity 

Rabc d = ftd a bc. d + (gca,Rb d — gcbRa, d — 5° aRbc + g d bRacj /2 

R 

+ (9ca9b d ~ gcb9a d ) ^ , (14) 

which defines the Weyl tensor. Expressing the Riemann tensor R a bc d in terms 
of the metric and its derivatives (or the Christoffcl quantities in a first order 
formalism) yields the desired equations. Note that for the physical Riemann 
tensor the vacuum Einstein equations imply R a bc d = C a bc d ■ 

We still miss differential equations for d a bc d and R a b- These can be obtained 
from the Bianchi identities ^7[ a Rb c ]d e , which in terms of the Weyl and tracefree 
Ricci tensors imply 

V d C abc d = (15) 
for the Weyl tensor of a vacuum spacetime (R a b = 0) and 

y b Ra b = \ Vai? • (16) 

While the Weyl tensor is conformally invariant, 
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this invariance does not hold for (|15[), instead however one can show that 

VdC a b c d = ft^d {d abc d ) , 

which implies 

V e d abc e = (17) 

if the vacuum Einstein equations hold in the physical spacetime. 
The Bianchi identity combined with the splitting ( [lij ) implies 

V a i? 6c - \7 b R ac = -— ((V a i?) g bc - (V 6 i?) g ac ) - 2 (V«jfl) d abc d • (18) 



The equations (11, I^ , |l^ , |r4| , |l7| , |lq ) then constitute the conformal field equa- 
tions for vacuum general relativity. Here the Ricci scalar R of g ab is considered 
a given function of the coordinates. For any solution (g a b,Rab,d abc d , ft), Rab is 
the traceless part of the Ricci tensor, and fi d a bc d the Weyl tensor of g ab . Note 
that the equations are regular even for ft = 0. 

The 3+1 decomposition of the conformal geometry can be carried out as 
usual in general relativity e.g. 

9ab = Kb - n a n b = ft 2 {Kb - n a h b ) , 

where h ab and h ab are the Riemannian 3-metrics induced by g ab respectively g ab 
on a spacelike hypersurface with unit normals n a , and equivalently n a = Q h a 
(our signature is (—>+,+>+))■ The relation of the extrinsic curvatures (k ab = 
\CnKb, k a b = \C n hab) is then easily derived as k ab = ft(k ab + ftoKh), where 
fta = n a V a ft. 

The additional variables R ab and d^ hc can be decomposed into spatial ob- 
jects by < 01 »i? a = n b h a c R bc , < 01 'i? ab = h a c h b d R bd , E ab = d e}cd h e a nf h c b n d , 
B a b = d* e j cd h e a n f h c b n d , where E ab and B ab are called the electric and mag- 
netic components of the rescaled Weyl tensor d abc d ■ 

Note that for regular components of h ab and k ab , the corresponding compo- 
nents of h ab and k ab with respect to the same coordinate system will in general 
diverge due to the compactification effect. However for the coordinate indepen- 
dent traces k = h ab k ab , k = h ab k ab of the extrinsic curvatures we get 

Qk = (k + 3ft ) , 

which can be assumed regular everywhere. Note that at J?, k = — 3J?o- Since 
is an ingoing null surface (with (V a i7)(V°^7) = but V a i7 ^ ), we have 
that fto < at J r+ . It follows that k > at ^ + . We will thus call regular 
spacelike hypersurfaces in M. hyperboloidal hypersurfaces, since in M. they are 
analogous to the standard hyperboloids t 2 — x 2 — y 2 — z 2 — p- in Minkowski 
space, which provide the standard example. Since such hypersurfaces cross 
but are everywhere spacelike in M. , they allow to access J? and radiation quan- 
tities defined there by solving a Cauchy problem (in contrast to a characteristic 
initial value problem which utilizes a null surface slicing). Note that in a glob- 
ally hyperbolic physical spacetime, hyperboloidal hypersurfaces will determine 
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the future of the physical spacetime, but not all of its past, we therefore call our 
studies semiglobal. 

The timelike vector t a = (d/dt) a is decomposed in the standard way into a 
normal and a tangential component: 

t a = Nn a + N a , N a n a = 0. (19) 

TV is called the lapse function, because it determines how fast the time evolution 
is pushed forward in the direction normal to S, and thus determines "how fast 
time elapses" . The tangential component N a , N a n a = 0, shifts spatial coordinate 
points with time evolution, accordingly N a is called shift vector. The lapse N 
and shift N a are not dynamical quantities, they can be specified freely and 
correspond to the arbitrary choice of coordinates: the lapse determines the slicing 
of spacetime, the choice of shift vector determines the spatial coordinates. 

We will not discuss the full 3 + 1 equations here for brevity, but rather 
refer to ||. Their most essential feature is that they split into constraints plus 
symmetric hyperbolic evolution equations ||. The evolution variables are h a b, 
k a b, the connection coefficients 7 a b c , <0,1) -Ra, <0,1) -Ra6, E a b, B a b, as well as J7, Qq, 
V a J?, V a V a /2 - in total this makes 57 quantities. In addition the gauge source 
functions q, R and N a have to be specified, in order to guarantee symmetric 
hyperbolicity they are given as functions of the coordinates. Here q determines 
the lapse as N — e q V det h and iV a is the shift vector. The Ricci scalar R can be 
thought of as implicitly steering the conformal factor fl. 

The constraints of the conformal field equations (see (14) of ||) are regular 
equations on the whole conformal spacetime {M. , g a b) , but they have not yet 
been cast into a standard type of PDE system, such as a system of elliptic PDEs 
(recently however, some progress in this direction has been achieved by Butscher 
[p0[). Therefore some remarks on how to proceed in this situation are in order. 
A possible resolution is to resort to a 3-step method [|np|,p0||: 

1. Obtain data for the Einstein equations: the first and second fundamental 
forms h a b and k a b induced on S by g a b, corresponding in the compactified 
picture to h a b, k a b and Q and Qq. This yields so-called "minimal data". 

2. Complete the minimal data on E to data for all variables using the conformal 
constraints - in principle this is mere algebra and differentiation. 

3. Extend the data from £ to S in some ad hoc but sufficiently smooth and 
"well-behaved" way. 

In order to simplify the first step, numerical implementations [|6|J^,|20|| so far 
have been restricted to a subclass of hyperboloidal slices where initially k a b is 
pure trace, k a b = \h a bk. The momentum constraint 

V b k ab - V„fe = Q (20) 

then implies k = const. ^ 0. We always set k > 0. In order to reduce the 
Hamiltonian constraint 

(3) i? + k 2 = ~k ab k ab 
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to one elliptic equation of second order, we use a modified Lichnerowicz ansatz 

h ab = fi~ 2 (jfhab 

with two conformal factors J? and <p. The principal idea is to choose h a b and i?, 
and solve for <j>, as we will describe now. First, the "boundary defining" function 
fl is chosen to vanish on a 2-surface S - the boundary of £ and initial cut of J? 
- with non- vanishing gradient on S. The topology of iS is chosen as spherical for 
asymptotically Minkowski spacctimes. Then we choose h a b to be a Ricmannian 
metric on U, with the only restriction that the extrinsic 2-curvature induced by 
hab on S is pure trace, which is required as a smoothness condition pif . With 
this ansatz h a b is singular at S, indicating that S represents an infinity. The 
Hamiltonian constraint then reduces to the Yamabe equation for the conformal 
factor (j): 

Af2 2 A(j>- 4n(y a n)(y a (t)) - Q (3) i?i? 2 + - 3(v Q j?)(V Q i?)^ = ip^ 5 . 

This is a semilinear elliptic equation - except at 5, where the principal part 
vanishes for a regular solution. This however determines the boundary values as 

4> 4 = 2- 2 {V a Q){V a Q) . (21) 

Existence and uniqueness of a positive solution to the Yamabe equation and the 
corresponding existence and uniqueness of regular data for the conformal field 
equations using the approach outlined above (assuming the "pure trace smooth- 
ness condition) have been proven by Andersson, Chrusciel and Friedrich [|T] . 

If the Yamabe equation is solved numerically, the boundary has to be chosen 
at S, the initial cut of J' , with boundary values satisfying (|2l|). If the equation 
would be solved on a larger grid (conveniently chosen to be Cartesian) , boundary 
conditions would have to be invented, which generically would cause the solu- 
tion to lack sufficient differentiability at S, see Hiibner's discussion in ||. This 
problem is due to the degeneracy of the Yamabe equation at S. Unfortunately, 
this means that we have to solve an elliptic problem with spherical boundary. 

The constraints needed to complete minimal initial data to data for all evolu- 
tion variables split into two groups: those that require divisions by the conformal 
factor J? to solve for the unknown variable, and those which do not. The latter do 
not cause any problems and can be solved without taking special care at fl — 0. 
The first group, needed to compute {1 - 1} R, E a b and B a b, however does require 
special numerical techniques to carry out the division, and furthermore it is not 
known whether solving them on the whole Cartesian time evolution grid actu- 
ally allows solutions which are sufficiently smooth across J? . Thus, at least for 
these we have to find some ad-hoc extension. There are however also examples of 
analytically known initial data, e.g. for the Minkowski and Kruskal spacetimes, 
where all constraints are solved on the whole Cartesian time evolution grid. 
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3 Examples: Different Ways to Compactify Minkowski 
Spacetime 

The examples presented in this section help to illustrate the compactification 
procedure - in particular its inherent gauge freedom, and they yield interesting 
numerical tests, some of which will be presented in Sec. pi 



3.1 Almost Static Compactification of Minkowski Spacetime 

From the perspective of hyperboloidal initial data, the simplest way to compact- 
ify Minkowski spacetime is to choose the initial conformal three-metric as the 
flat metric, h a b — 5 a b, to set k a b = h a b, which solves the momentum constraint 
( ^0|) , and to choose the conformal curvature scalar R g ^ as spherically symmet- 
ric, R g = R g (x 2 + y 2 + z 2 ). We know from (3l) that a unique solution to the 
constraints exists, it is not hard to see that it has to be spherically symmetric. 
Furthermore, it is topologically trivial. From Birkhoff's theorem we can thus 
conclude that we are dealing with Minkowski spacetime. Choosing the simplest 
gauge q — 0, N a — 0, R g — 0, the resulting unphysical spacetime is actually 
Minkowski spacetime in standard coordinates: 

ds 2 = -At + AS 2 = tt 2 (-AT 2 + AR 2 + R 2 {AO 2 + sin 2 6>d</> 2 ) , 

where dS 2 is the standard metric on R 3 , dS 2 — dr 2 + r 2 (dO 2 + sin 2 #d</> 2 ), and 
the conformal factor is 

f}=(R 2 -T 2 y 1 = {r 2 -t 2 ) , (22) 

where 

R T 

r ~R 2 -T 2 ' *~ R 2 -T 2 ' 

This setup has been chosen as the basis of Hubner's numerical study of weak 
data evolutions J^] . With the initial cut of y + chosen at x 2 + y 2 + z 2 = 1 , i + 
is located at coordinate time t = 1, the generators of being straight lines at 
an angle of 45°. 

This conformal representation of Minkowski spacetime constitutes an "al- 
most static" gauge - since the spatial geometry is time- independent, so are all 
evolution variables except for the conformal factor Q. The physical region inside 
of contracts to the regular point i + within finite time. This feature is shared 
with the standard "textbook" example of conformally compactifying Minkowski 
spacetime, which takes the form of a map into part of the Einstein static universe 
with R g = 6, 

ds 2 = -At 2 + AS 2 = Q 2 (-AT 2 + AR 2 + R 2 (AO 2 + sin 2 <9d</> 2 )) , (23) 

3 We change notation from R to R g for this section to avoid confusion with a coordinate 
R we will introduce below. 
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where d£ 2 is the standard metric on S* 3 , dS 2 = dg 2 + sin 2 g (d8 2 + sin 2 9d(j) 2 ), 
and the conformal factor is 

n 2 = 4 (1 + (T - R) 2 )- 1 (1 + (T + i?) 2 )- 1 = 4cos 2 ^ cos 2 ^ . 

Here the coordinate transformations are 

£> = arctan(T + i?) - arctan(T - R) , (24) 
f = arctan(T + R) + arctan(T - R) . (25) 

In these coordinates Minkowski spacctime corresponds to the coordinate ranges 

-7T <t + g<ir , (26) 
-7T <t-g<ir, (27) 
e > o ■ (28) 

For details and pictures of this mapping see the discussions by [Q, || or |17| . 
Alternatively, we can choose stereographic spatial coordinates such that 

dS 2 = uj 2 (dr 2 + r 2 (d6 2 + sin 2 9dcj) 2 )) , u - 2/(1 + r 2 ) . 

or we may absorb the spatial conformal factor into the spacetime conformal 
factor by rescaling to 

ds' 2 = -uj- 2 dt 2 + dr 2 + r 2 (d<9 2 + sin 2 dd<p 2 ) , (29) 

which yields the lapse to be N = 1 (q = — 31ogw), respectively N = uj~ 2 
(q = — 2 logo;). Note that in the numerical code we use Cartesian coordinates 
x = r sin 9 cos 4>, y — r sin 9 sin <p, z — r cos 9. 

The conformal transformation leading to (^9|) changes the scalar curvature 
from R g — 6 to R g = —12 (1 + r 2 )~ 2 . We will see below in Sec. || that these 
simple variations in gauge source functions and conformal rescaling lead to nu- 
merical representations which are quite different, e.g. with regard to accuracy 
and robustness. 



3.2 A Static Hyperboloidal Gauge for Minkowski Spacetime 

By translating a standard hyperboloid in Minkowski spacetime along the tra- 
jectories of the d/dt Killing vector, one can obtain a gauge where not only the 
conformal spacetime is static, but also the conformal factor is time-independent 
- thus also the physical geometry and all evolution variables of the conformal 
field equations can be made time independent (this has been pointed out me by 
M. Weaver, and I essentially follow her notes below). See also a talk given by 
V. Moncrief jig] , that we have become aware of after starting to work with this 
gauge. 

In this gauge the point i + is not brought into a finite distance and remains in 
the infinite future. This conformal gauge is particularly useful for stability tests. 
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To derive this static metric we start with spherical coordinates (T, R, 9, </>) 
on Minkowski space, where the metric is 

dS 2 = -dT 2 + di? 2 + R 2 (dd 2 + sin 2 9d(f> 2 ) . (30) 

A family of standard hyperboloids with time translation parameter t is given by 

(T - tf - R 2 = 1 . 

We transform now to new coordinates (t,g,9,</>), where the level surfaces of t 
are the standard hyperboloids, and g(R) is chosen as a new radial parameter 
on the hyperboloids. Setting T — t + cosh g and R = sinh g, the physical metric 
becomes 

dS 2 = -d< 2 - 2 sinh g dg dt + dg 2 + sinh 2 g (d9 2 + sin 2 6d<j) 2 ) . (31) 

For simplicity we choose the conformal three-metric to be flat and introduce new 
spherical coordinates (r,9,<fi), such that 

ds 2 \ t=const . = dr 2 + r 2 (d9 2 + sin 2 9dtf) . (32) 
Since — f2 2 h a b we get Q = dr/dg and 

. ^ = /'^ 

The limits of integration are given by the fact that lim^oo r — 1. Performing 
the integrals one finds that 



e e - 1 



l^l + R 2 -!) (34) 



ee + 1 R 
and 

1-r 2 

Q = —Ji- , (35) 

our choice thus maps to the timelike cylinder r = 1. 

The computer time coordinate t is a Bondi time coordinate on ^ . In coor- 
dinates (t, r, 6, (f>) the conformal metric reads 

„2 



ol- 



ds 2 = -f2 2 dt 2 - 2rdrdt + dr 2 + r 2 (d6» 2 + sin 2 9d(j) 2 ) , (36) 



ds 2 = -n 2 dt 2 -2dt (xdx + ydy + zdz) + dx 2 + dy 2 + dz 2 (37) 



in Cartesian coordinates x — rsin0cos0, y — rsin^sin^, z — rcos9, which are 
used in the numerical code. The shift vector is thus given by N l = —x 1 , and 
the lapse can be computed from — N 2 + h a bN a N b = g u (as implied by ([l9|)) as 
iV = (1 + r 2 )/2. The three metric has unit determinant, so q = InN. Note that 
the shift vector does not become "super luminal" beyond because the lapse 
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is growing faster than the shift - gtt is nonnegative everywhere, and zero only 
at J? + . The conformal Ricci scalar is 



_ (i-r2)( 3 + r») 

- 12 (l + r 2 ) 3 ' ( ^ 



which vanishes at 

For a numerical calculation one needs the minimal initial data set, 

(h ab ,n,k ab ,O ) (39) 

and the gauge source functions, (R, N, N a ). In a numerical calculation in which 
the Yamabe equation is solved to find fi, one gives (h a b, 0,trk). In a test case 
such as this which is an explicitly known solution, one can just take Q = Q. It 
remains therefore to calculate k a b and fl$. From 

k a b = -^{dthab - C N h ab ) (40) 

we find that the components of the extrinsic curvature are ku = -h^ij ancl 



y — n y 

k = jj. From the identity k — Q k — 3 £2® we find that 



2r 2 

^ = ■ (41) 



4 History 

This section tries to give a broad overview of what has been achieved so far in the 
field of numerical treatment of the conformal field equations. Historically, this 
field was started by Peter Hiibner by studying a scalar field coupled to gravity in 
spherical symmetry in his PhD thesis ]3^| finished in 1993. His subsequent work 
has lead to the development of both a 2D and a 3D evolution code, formulated 
in "metric" variables. Jorg Frauendiener has also developed an independent 2D 
code, formulated in frame variables. 



4.1 Early work in spherical symmetry 

The first numerical implementation of the conformal field equations is due to Pe- 
ter Hiibner, who has studied the spherically symmetric collapse of scalar fields in 
his PhD thesis ||^] , and subsequently in |}3| . In his gauge both future null infinity 
and future timelike infinity are compactified, and the whole spacetime 
is covered in finite coordinate time. Hiibner studies the global structure of the 
spacetime, including the appearance of singularities and the localization of the 
event horizon. To handle the latter, floating point exceptions are caught and grid 
points are flagged as "singular" , grid points whose values depend on information 
from singular grid points are correspondingly flagged as singular as well. Even 
though this method does not allow to actually trace the singularity in a strict 
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sense (computers can not actually deal with infinite values), the method traces 
the singularity as tightly as possible. In contrast to typical black hole excision 
schemes, which are based on locating the apparent horizon, this scheme could 
thus be termed "tight excision" . The method has not yet been implemented in 
higher dimensions, where one has to face more intricate technical problems, and 
where the structure of the singularity is likely to be much more complex as well. 
The paper also studies critical collapse. Hiibner's results are consistent with the 
black hole mass power-law scaling with the correct exponent, however no echoing 
related to discrete self-similarity has been seen in his results. This has created 
some discussion, whether the results of other authors are numerical artifacts, 
or artefact's of boundary conditions at finite distance, however numerical criti- 
cal collapse simulations in a compactified characteristic framework have recently 
shown both the correct power-law scaling and discretely self-similar echoing Q . 

The coordinates in this approach are based on the geometric structure of dou- 
ble null-coordinates available in spherical symmetry. Unfortunately this choice 
does not generalize in the absence of spherical symmetry. Finding a gauge that 
would allow to run, say, the Kruskal spacetime in a 3D code for "arbitrarily 
long" Bondi times, is an open problem, where significant insight could be gained 
from studying more general gauges in a manifestly spherically symmetric code. 



4.2 Axially symmetric spacetimes with toroidal Scri in the frame 
formulation 

Following Hiibner's encouraging results for spherically symmetric simulations 
[[32j|3j|, numerical codes have been developed by Frauendiener and Hiibner to 
study axially symmetric spacetimes. For simplicity, e.g. to avoid numerical sta- 
bility problems at the axis of symmetry, and to avoid problems associated with a 
J* of spherical topology - which does not align with Cartesian coordinates - both 
Hiibner and Frauendiener considered the asymptotically A3-spacetimes p5|p^ |, 
which do not possess an axis of symmetry and where J 1 has toroidal topology. 
These spacetimes are modeled after the A3-metric in the Ehlers-Kundt classi- 
fication [0, which provides an analogue of the Schwarzschild metric in plane 
symmetry. These spacetimes are not physical, but they contain a large class 
of nontrivial radiative vacuum spacetimes, which make them an interesting toy 
model to study numerical techniques, gauges, and the extraction of radiation Q 
These axisymmetric codes thus have been designed to treat the vacuum case, and 
matter couplings have not yet been implemented. An advantage for code-testing 
is that exact solutions are known |35|j4l|| . 

In the first |l^] of a series of papers (l^Jl^,^,^l) on his axisymmetric code, 
Frauendiener gives a nice overview of the motivation for using the conformal 
field equations of numerical simulations of isolated systems. He discusses the 
conformal field equations in the space spinor formalism pis] , which is chosen 

4 One of the notable differences to the Minkowski case is that one can only define a 
Bondi-mass but no Bondi four-momentum. 



18 Sascha Husa 



because of compactness of notation, and because it allows a very straightfor- 
ward 3+1 decomposition of the equations, rendering the equations in symmetric 
hyperbolic form. His formalism contains 8 free functions which determine the 
gauge: the harmonicity F :— V ' C VH determines the choice of time coordinate t, 
the shift is given in terms of frame coefficients, the scalar curvature R (A in his 
notation) of the compactified spacetime, and an imaginary and symmetric space 
spinor field Fab (he. three numbers) which determines rotations of the spatial 
frame (for Fab = the frame is transported via Fermi- Walker transport). He 
also discusses the implications of the assumptions of the toroidal symmetry, in 
particular for the choice of gauge - e.g. the adaption of the frame. 

In the second paper |19|| of the series Frauendiener discusses his numeri- 
cal methods and gauge choices, and presents results from evolutions of initial 
data corresponding to the exact solution presented in |35| . Here one of the two 
Killing vectors is disguised by a coordinate transformation. The numerical evo- 
lution proceeds via a generalization of the Lax-Wendroff scheme to 2D, which 
Frauendiener proves to be stable and second order accurate. The time step is 
such that the numerical domain of dependence is contained in the domain of 
dependence as defined by the equations. An essential difficulty - as usual - 
is posed by the treatment of the boundary. Well-posedness of the associated 
initial-boundary-value problem has not yet been proven, and numerical analysis 
can only provide rough guidelines to work out stable algorithms |3^] . Frauendi- 
ener's boundary treatment is based on the identification of ingoing and outgoing 
modes at the boundary, as determined from the symmetric hyperbolic character 
of the equations. He sets boundary values for inward-propagating quantities (e.g. 
motivated by the exact solution), and sets values for the outward propagating 
quantities by extrapolation from the interior. This method can be applied just 
a few grid-points outside of and is found to be stable as long as the gauge 
source functions dot not depend on the evolution variables - which would change 
the characteristics. Note that the constraints will in general not be satisfied on 
the boundary, which may trigger constraint-violating modes of the equations. 

Frauendiener gives a detailed discussion of the problems associated with the 
choice of the gauge, and performs a number of numerical experiments in this 
respect, evolving data corresponding to exact solutions J3||,0 with singular i + . 
One of the problems is, that if the gauge source functions are allowed to depend 
on the evolution variables, this will change the characteristics of the system 
and will in general spoil the symmetric hyperbolic character of the system. Ex- 
periments in this direction, where F = F(N, K) indeed exhibited a boundary 
instability. Regarding the choice of time coordinate, that is the harmonicity func- 
tion F, several choices are tested: a "natural" gauge, which is taken from the 
exact solution, the Gauss gauge (where the lapse N is spatially constant), the 
harmonic gauge, F = 0, and a family of gauges that interpolates between the 
"natural" and harmonic gauge. 

The "natural" gauge is found to provide the best performance and the ap- 
proach to the singularity is found to be essentially limited by machine precision. 
The harmonic gauge leads to a coordinate singularity before reaching the sin- 
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gularity, this feature is shared by most of the gauges that interpolate between 
natural and harmonic gauge. For the "Gauss" gauge with N = const, caustics 
(coordinate shocks) develop quickly and crash the simulations. 

Regarding the choice of shift vector, a prescription for J" fixing - that is 
steering the evolution of the surface Q = is discussed, which can be easily 
implemented in Frauendiener's formulation. This however relies on the specific 
form of the frame equations, and does not carry over to to equations as used in 
Hubner's codes [||. In particular he studies the case of freezing" - holding 
the coordinate position of in place such that no loss of resolution occurs in 
the physical domain. 

Finally, he discusses the extraction of gravitational radiation, e.g. by com- 
puting the Bondi mass, and shows some results. 

In order to study more general spacetimes, Frauendiener has implemented 
a numerical scheme for determining hyperboloidal initial data sets for the con- 
formal field equations by using pseudo-spectral methods as described in 
He uses the implicit approach of first solving the Yamabe equation, and then 
carrying out the division by the conformal factor for certain fields which vanish 
on J' ' . The challenge there is to numerically obtain a smooth quotient. The di- 
vision problem is treated by a transformation to the coefficient space, where a 
QR-factorization of a suitable matrix is used, and then transforming back. 

In H^] Frauendiener gives a pedagogical discussion of the issue of radiation 
extraction in asymptotically flat space-times within the framework of conformal 
methods for numerical relativity. The aim is to show that there exists a well de- 
fined and accurate extraction procedure which mimics the physical measurement 
process and operates entirely intrinsically within . The notion of a detector 
at infinity is defined by idealizing local observers in Minkowski space. A detailed 
discussion is presented for Maxwell fields and the generalization to linearized 
and full gravity is performed by way of the similar structure of the asymptotic 
fields. 

Recently, Hein has written a 2D axisymmetric code that allows for an axis 
[E2J, i.e. can treat the physical situation with a of spherical topology. The 
usual problem of the coordinate singularity at the axis in adapted coordinates 
was solved by using Cartesian coordinates, following a method developed by 
Alcubierre et al. The code has so far been tested by evolving Minkowski 
spacetime in various gauges, further tests with nontrivial spacetimes are cur- 
rently underway. 

4.3 Metric-based 2D and 3D codes 

The basic design of Hubner's approach is outlined in [^J, where he presents the 
first order time evolution equations as obtained from a 3+1 split of the con- 
formal field equations. The evolution equations can be brought into symmetric 
hyperbolic form by a change of variables. He discusses his motivation of avoid- 
ing artificial boundaries and how the conformal field equations formally allow 
placement of the grid boundaries outside the physical spacetime. 
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A particularly subtle part of the evolution usually is the boundary treatment. 
In the conformal approach we are in the situation that the boundary can actually 
be placed outside of the physical region of the grid - this is one of its essential 
advantages! In typical explicit time evolution algorithms, such as our Runge- 
Kutta method of lines, the numerical propagation speed is larger than the speed 
of all the characteristics (in our case the speed of light). Thus does not shield 
the physical region from the influence of the boundary - but this influence has 
to converge to zero with the convergence order of the algorithm - fourth order in 
our case. In principle one therefore does not have to choose a "physical" bound- 
ary condition, the only requirements are stability and "practicality" - e.g. the 
boundary condition should avoid, if possible, the development of large gradients 
in the unphysical region to reduce the numerical "spill over" into the physical 
region, or even code crashes. It seems likely however, that this practicality re- 
quirement will eventually lead to a treatment of the boundary which satisfies 
the constraints at the boundary. 

Hiibner develops the idea of modifying the equations near the grid boundaries 
to obtain a consistent and stable discretization. The current implementation of 
the boundary treatment relies on this introduction of a "transition layer" in the 
unphysical region, which is used to transform the rescaled Einstein equations 
to trivial evolution equations, which are stable with a trivial copy operation at 
outermost gridpoint as a boundary condition (see for details and references) . 
He thus replaces 

d t f + A l dJ -6 = 

by 

d t f + a(Q) (A l dJ -b)=0, 

where a is chosen as a(Q) = for fl < i?o < tt\ < and 1 for fi > Q\. 
One potential problem is that the region of large constraint violations outside 
of J? may trigger constraint violating modes of the equations that can grow 
exponentially. Another problem is that a "thin" transition zone causes large 
gradients in the coefficients of the equations - thus eventually leading to large 
gradients in the solution, while a "thick" transition zone means to loose many 
gridpoints. If no transition zone is used at all, and the Cartesian grid boundary 
touches J? ', the ratio of the number of grid points in the physical region versus 
the number of grid points in the physical region is already 7r/6 ~ 0.52. 

Furthermore he discusses his point of view concerning possible advantages 
of the conformal approach and discusses potential problems of the Cauchy and 
Cauchy- Characteristic matching approaches to numerical relativity. He outlines 
the geometric scenario of his approach and stresses that these techniques allow, 
in principle, to calculate the complete future of scenarios such as initial data for 
./V black holes. 

The second paper || of the series deals with the technical details of construc- 
tion of initial data and of the time-evolution of such data. The second and fourth 
order discretizations, which are used for the construction of the complete data 
set and for the numerical integration of the time evolution equations, are de- 
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scribed and their efficiencies compared. Results from tests for A3 and disguised 
Minkowski spacetimes confirm convergence for the 2D and 3D codes. 

The simplest approach to the division by J? would be an implementation of 
l'Hospital's rule, however this leads to nonsmooth errors and consequently to a 
loss of convergence [0 . Instead Hiibner [|| has developed a technique to replace 
a division g = f / f2 by solving an elliptic equation of the type 

V a V a (Q 2 g~ Qf) = 

for g (actually some additional terms added for technical reasons are omitted 
here for simplicity). For the boundary values f2 2 g — Qf = 0, the unique solution 
is g = f / fi. The resulting linear elliptic equations for g are solved by the same 
numerical techniques as the Yamabe equation. For technical details see Hiibner 

Finally, we have to extend the initial data to the full Cartesian spatial grid 
in some way. Since solving all constraints also outside of £ will in general not 
be possible in a sufficiently smooth way || , we have to find an ad hoc extension, 
which violates the constraints outside of J* but is sufficiently well behaved to 
serve as initial data. The resulting constraint violation is not necessarily harmful 
for the evolution, since J causally disconnects the physical region from the 
region of constraint violation. On the numerical level, errors from the constraint 
violating region will in general propagate into the physical region, but if our 
scheme is consistent, these errors have to converge to zero with the convergence 
order of the numerical scheme (fourth order in our case). There may of course 
still be practical problems that prevent us from reaching this aim: making the 
ad- hoc extension well behaved is actually quite difficult, the initial constraint 
violation may trigger constraint violating modes in the equations, which take 
us away from the true solution, singularities may form in the unphysical region, 
etc. 

Since the time evolution grid is Cartesian, its grid points will in general 
not coincide with the collocation points of the pseudo-spectral grid. Thus fast 
Fourier transformations can not be used for transformation to the time evolution 
grid. The current implementation instead uses standard discrete ("slow") Fourier 
transformations, which typically take up the major part of the computational 
effort of producing initial data. 

It turns out, that the combined procedure works reasonably well for certain 
data sets. For other data sets the division by fi is not yet solved in a satisfactory 
way, and constraint violations are of order unity for the highest available resolu- 
tions. In particular this concerns the constraint V b E a b = - { ^ a bck bd B d c ((14d) 
in H ) , since E a i> is computed last in the hierarchy of variables and requires two 
divisions by J?. Further research is required to analyze the problems and cither 
improve the current implementation or apply alternative algorithms. Ultimately, 
it seems desirable to change the algorithm of obtaining initial data to a method 
that solves the conformal constraints directly and therefore does not suffer from 
the current problems. This approach may of course introduce new problems like 
an elliptic system too large to be handled in practice. 
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The time evolution algorithm is an implementation of a standard fourth order 
method of lines (sec e.g. |f(|). In the method of lines we formally write 

d t f = B(f,dif), (42) 

where B(f, dif) = —A l (f)dif + b(f). Discretizing the spatial derivatives param- 
eterizes the ordinary differential equations by grid point index. For the present 
code, fourth order accurate centered spatial differences have been implemented, 
e.g. for the rr-derivative as: 

d x f — > \2Ax ^~f i+2 >j> k 8/;+i,j,fc — &ft-i,j,k + fi-i,j,k) ■ 

The numerical integration of the ordinary differential equations proceeds via 
the standard fourth order Runge-Kutta scheme: 

tl+i _ fl + lful 4- 2k l+1/i 4- 2k l+1/2 + k l+3/i \ (41) 

where 

K,j,k = AtB(f- j k ,dif- j k ) , 

,1 + 1/4 _ A , p, ,1+1/4, r\ A+l/is ,2+1/4 _ A , 1 , J 

J+l/2 _ A , „(. 1+1/2 p. f l + l/2, A+l/2 _ A 1 j. 1+1/4 

K iJ,fc —^ Aln \Ji, 3 ,k > iJi,j,k )' h,j,k — Ji,j,k 2 ' 

■ i+3/4 _ . R , ,J+3/4 fl ,i+3/4, ,i+3/4 _ jJ+1/2 

rc i,i,fe ^ ASr -°Ui J ,fc >°iJi,j,k /' ■'i.i.fe — Ji,j,k + K i,j,k ■ 

Additionally, a dissipation term of the type discussed in theorems 6.7.1 and 6.7.2 
of Gustafsson, Kreiss and Oliger jfO| is added to the right-hand-sides to damp out 
high frequency oscillations and keep the code numerically stable. The dissipation 
term used is 0Q2 := (Ax) 5 YliLi 9i 6 f, where the spatial derivatives are 
discretized as 

d x & fl 3 ,k -> r^)6 (fi-3,j,k " 6 /i-2,j,fe + 15 /'-lj,fc 
_ 20/- j fe + 15/- +1 j fe — 6/j_|_ 2 ,j,)fc + fi+3,j,k) ■ 

Numerical experiments show that usually small amounts of dissipation (ct of 
order unity or smaller) are sufficient, and do not change the results in any signif- 
icant manner. Numerical tests for Minkowski spacetime with disguised symme- 
tries and an explicitly known A3-like solution with radiation p4| are described 
in §. 

Further extensive tests of the 2D code have been performed by Weaver [Q . 
She studied the choice of gauge source functions for an A3-like solution, solving 
the Yamabe equation for the conformal factor. She found that for this solu- 
tion it is quite simple to prescribe a shift so that .y is fixed to a very good 
approximation. She also studied use of the gauge source function q to prolong 
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the numerical simulation inside physical spacetime. In cases where q = re- 
sults in a "singularity" developing outside physical spacetime (which causes the 
code to crash), prescription of q so that the evolution inside physical spacetime 
is prolonged compared to outside allows the simulation to essentially cover the 
physical spacetime to the future of the initial data surface. She thus found that 
in this context the ad hoc prescription of gauge source functions was sufficient 
to achieve desired effects, and caused no instabilities. Also she explored the ef- 
fect of turning off the transition zone, while still simply copying data at the 
outer grid boundary into the ghost zone, along with prescription of q so that 
the evolution is slowed down at the outer boundary. In the A3-like 2D runs this 
alternative boundary treatment was successful, and avoided problems created 
by the transition zone. 

In the third part of the series (tJ , a pseudospectral solver for the constraints is 
described. Since the implementation depends on the topology, it discusses both 
the asymptotically A3 and asymptotically Minkowski cases. At the end also some 
remarks are made about a possible extension to the multi-black-hole case, using 
a multi-patch scheme (the Schwarz alternating procedure). 

In the fourth part of the series || Hiibner presents results of 3D calculations 
for initial data which evolve into a regular point i + , and which thus could be 
called "weak data" . The initial conformal metric is chosen in Cartesian coordi- 
nates as 

ds 2 = ( 1 + (x 2 + 2y 2 )^j dx 2 + dy 2 + dz 2 . (44) 

The boundary defining function fl appearing in this ansatz is chosen as fl = 
(l — {x 2 + y 2 + z 2 )) /2, it is used to satisfy the smoothness condition for the 
conformal metric at J?. For the gauge source functions Hiibner has made the 
"trivial" choice: R — 0, N a = 0, q = 0, i.e. the conformal spacetime has vanishing 
scalar curvature, the shift vanishes and the lapse is given by N — e q Vdeth = 
Vdet h. This simplest choice of gauge is completely sufficient for A = 1 data, and 
has lead to a milestone result of the conformal approach - the evolution of weak 
data which evolve into a regular point i + of Ai , which is resolved as a single grid 
cell. With this result Hiibner has illustrated a theorem by Friedrich, who has 
shown that for sufficiently weak initial data there exists a regular point i + of M. 
[ pi) . The complete future of (the physical part of) the initial slice can thus be 
reconstructed in a finite number of computational time steps. This calculation is 
an example of a situation for which the usage of the conformal field equations is 
ideally suited: main difficulties of the problem are directly addressed and solved 
by using the conformal field equations. 

A natural next question to ask is: what happens if one increases the amplitude 
Al To answer this question, I have performed and analyzed runs for integer 
values of A up to A = 20, preliminary results have been presented in |47|. While 
for A — 1,2 the code was found to be able to continue beyond i + without 
problems, for all higher amplitudes the "trivial" gauge leads to code crashes 
before reaching i + . While the physical data still decay quickly in time, a sharp 
peak of the lapse develops outside of J" and crashes the code after Bondi time 
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- 8(320M) for A = 3 and - 1.5(3M) for A = 20 (here M is the initial Bondi 
mass). A partial cure of the problem was obtained using a modified gauge source 
function q = —r 2 /a (N = e~ r '°\/det h), where a is tuned such that one gets 
a smooth lapse and smooth metric components. For, A — 5, for example, a 
value of a = 1 was found by moderate tuning of a (significantly decreasing or 
increasing a crashes the code before the regular i + is reached). Unfortunately, 
this modification of the lapse is not sufficient to achieve much higher amplitudes. 
As A is increased, the parameter a requires more fine tuning, which was only 
achieved for A < 8. For higher amplitudes the code crashes with significant 
differences in the maximal and minimal Bondi time achieved, while the radiation 
still decays very rapidly. Furthermore, the curvature quantities do not show 
excessive growth - it is thus natural to assume that we are still in the weak- 
field regime, and the crash is not connected to the formation of an apparent 
horizon or singularity. While some improvement is obviously possible through 
simple non-trivial models for the lapse (or other gauge source functions), this 
approach seems quite limited and more understanding will be necessary to find 
practicable gauges. An interesting line of research would be to follow the lines 
of [ff!)| in order to find evolution equations for the gauge source functions which 
avoid the development of pathologies. 

Schmidt has presented hyperboloidal initial data for the Kruskal spacetimc, 
a hyperboloidal foliation for the future of these hyperboloidal initial data |Q 
and results from numerical simulations evolving these initial data with differ- 
ent gauges, which have been performed by Weaver with Hubner's 3D code. The 
explicit hyperboloidal version of the Kruskal spacetime is very useful for numer- 
ically testing the conformal approach in the treatment of black hole spacetimes. 
These runs have been performed in octant mode. The runs typically proceed 
until the determinant of the three metric becomes negative H] , caused by some 
feature in the exact solution which is no longer adequately resolved and which is 
growing, leading to large narrow spikes in the numerical data. Future work will 
have to be directed toward improving the choice of gauge source functions such 
that rapidly growing sharp features are avoided. 

In the next section, I will present new results obtained with the 3D code for 
asymptotically Minkowski spacetimes, which will illustrate some of the current 
problems. One of these is the presence of exponentially growing constraint vio- 
lating modes. The problem of controlling the growth of the constraints for the 
conformal field equations has first been addressed by Florian Siebel in a diploma 
thesis |5Cj , and subsequently by Hiibner and Siebel in |Q . The key idea in this 
work is to develop a A-system |35| for the conformal field equations in 1+1 
dimensions (with toroidal ^'s). A A-system is an enlarged evolution system, 
where evolution equations for the constraints are added in, consistently with 
symmetric hyperbolicity. One then has a large parameter space of coefficient 
functions available, in which to find choices such that the new system has the 
constraint surface as an attractor. The main conclusion of this work is that it was 
not possible to significantly improve the fidelity of the numerical calculations. 
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In those cases where moderate improvements regarding the constraints could be 
achieved, the deviation from the known exact solution would get larger. 

5 Results from 3D calculations 

All the results presented in this section have been performed with 121 3 grids on 
32 processors of the AEFs SGI origin 2000. The outer boundary has been placed 
at a radius of r = 1.15 in these runs is initially located at a radius r — 1). 



5.1 Minkowski data 

We will first discuss some results for Minkowski spacetime, which in spite of its 
simplicity provides some nontrivial numerical tests. As has been first demon- 
strated by Hubner in || , for weak data - in particular Minkowski space - it is 
possible via the conformal approach to cover the whole domain of dependence 
of initial data reaching out to with a finite number of time steps. Let us 



thus first consider the gauges of Sec. 3.1, where the compactified geometry is 



time-independent, but a time-dependent conformal factor Q is responsible for 
contracting the cuts of ^ + to a point within finite coordinate time. 

We have compared the gauges where the conformal spacetime is Minkowski, 
22]), the Einstein static universe (p3|), or the spacetime given by (p^). Essentially, 
the result is that the Minkowski case yields the highest accuracy, the Einstein 
universe case works in principle, and in the case ( p9| ) the code crashes before 
reaching i + . In Fig. |l| the Minkowski and Einstein universe cases are compared 




Fig. 1. Comparing the Minkowski (solid line) and Einstein universe (dashed line) cases: 
left the value of the constraint V^J? = Q x at the center is plotted versus coordinate 
time, in the right image h xx — 1 is plotted vs. coordinate time (where t is scaled such 
that t(i + ) = 1) 



by plotting h xx — 1 and the value of the constraint \7 x f2 — fl x at the center 
versus coordinate time (where t is scaled such that t(i + ) = 1. The Minkowski 
case - denoted by the unbroken line - clearly yields better accuracy, although 
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the growth of h xx — 1 is faster, and approximately exponential during the later 
stage of "physical" evolution. Note that the constraint grows very fast in both 
cases. 

Figures || - J5 show a comparison of the less optimal Einstein universe case 
with the case (E9) to illustrate some of the problems one expects in the evolution 
of nontrivial spacetimes. Fig. ^ shows the time evolution of h xx along the positive 
x-axis versus coordinate time for the Einstein universe case and for the case of 
(|29|). Fig. H compares the corresponding contour lines. While no deviation from 
staticity is visible for the Einstein universe case, the other case shows a rapidly 
growing peak in h xx and the lapse (shown in Fig. ||) (and thus of det h) , which 
is located in the transition zone outside of J r+ . Eventually this feature can not 
be resolved any more, and the code crashes. In the Einstein static case the code 
was simply stopped by running out of time in the queue. Fig. || shows the sum 
over the L 2 -norms (taken in the physical region) of all the constraints versus 
time. While in the Einstein static case the constraints show a rapid decrease in 
the physical region, followed by a steep growth after passing through i + , the 



case (E9[) exhibits roughly exponential overall growth almost from the start. 





Fig. 2. The value of the metric component h xx for x > is plotted versus coordinate 
time, the left image shows the Einstein universe case, the right image shows case (p9|), 
there the maximum of h xx in the region where S7 > is approximately at the value 5 



Results for the completely static gauge given by j37| ) are shown in Figs. || 
- ||. This gauge poses a harder challenge than the previous ones, where i + is 
reached in finite time. Now the goal is to maintain an indefinite stable evolution. 
However, the evolution shows exponential growth, illustrated in Figs. [7] and 
^| by the values of h xx and constraints V x h xx and V x f2 — fl x . Interestingly, 
however, the curvature invariants / and J are decreasing during the evolution 
as shown in Fig. [(]. The exponential blowup crashes the code at t ~ 5.1, this 
time seems to be roughly independent of resolution, size of time step, amount 
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Fig. 3. Contour lines of the metric component h xx for x > are plotted versus coordi- 
nate time, the left image shows the Einstein universe case, the right image shows case 
(||). The thicker line marks Q = 0, i.e. 




Fig. 4. The value of the lapse N for x > is plotted for the case (|29|) versus coordinate 
time, the left image shows the points where Q > 0, the right image shows all points 




Fig. 5. The sum over the L 2 -norms (taken in the physical region) of all the constraints 
is plotted versus coordinate time for the Einstein universe (left) and case ( p9| ) (right) 
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of dissipation, location of the boundary and location of the transition zone. A 
possible explanation are exponentially growing constraint violating modes on 
the continuum level. 




Fig. 6. The real parts of the curvature invariants I (left) and J (right) are plotted 
versus coordinate time for the static gauge of (|37]). The solid line is for the gridpoint 
at the center of the grid, the dashed line for a grid point at x = 0.996, y = z — 0, 
multiplied by a factor of 10 -5 for I and 10 -8 for J 




Fig. 7. The value of the metric component h xx for x > is plotted versus coordinate 
time with linear (left) and logarithmic (right) scaling for the static gauge of (^). 
Approximately exponential growth is obvious, the largest amplitude of the growth is 
in the center 



5.2 "Brill" data 

We use an axisymmetric Brill-wave type ansatz to look at initial data that 
contain radiation and set 

ds 2 = oj 2 (e 2Q (dg 2 +dz 2 ) + g 2 d(f 2 ) , (45) 
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where g 2 = x 2 + y 2 . With Q — \ ln(l + Af2 2 g 2 / (g 2 )) , in Cartesian coordinates 
the conformal three-metric becomes 

/ '1 + Ax 2 _Q 2 f AxyQ 2 f 
h B =uj 2 \ AxyQ 2 f Ay 2 f2 2 f 0_ 

\ l + AQ 2 f 

The axial symmetry makes it easier to analyze the data and choose the gauges. 
Here we set u = f = 1 and A = 1. 

Fig. shows the real part of the physical curvature invariant / = J? 6 / and 
the mass loss Mb- The curvature invariant I is computed both as a perturbation 
of the Einstein universe and case (p9|) (triangles) for a "Brill wave" with A = 1, 
to demonstrate the the physical initial data are indeed identical. The mass loss 
Mb is computed as a perturbation of the Einstein static case (R g = 6) and 
plotted in a logarithmic scale. Note that the falloff levels off at late times to a 
constant value due to numerical error. Note also that oscillations, like present 
here, are absent from the initial data corresponding to (^J) as shown in Fig. 5 

of gg. 

6 Conclusions and Outlook 

Bringing the conformal approach to numerical relativity to full fruition such 
that it can be used as a tool to explore new physics - in particular in black hole 
spacetimes - will be a long term effort. In order to contemplate the scope of this 
project, let us give a drastically oversimplified definition of the art of numerical 
relativity as a procedural recipe: 
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Fig. 9. The left image shows the real part of the physical curvature invariant I = fi ' I, 
computed as a perturbation of the Einstein universe (line) and case ((is]) (triangles) for 
a "Brill wave" with A = 1. The right image shows the corresponding mass loss function 
Mb, computed as as a perturbation of the Einstein static case (R g — 6) 



1. Find a well posed formulation of the initial(-boundary) value and initial data 
(constraint) problems for general relativity (optimally, well-posedness should 
be a theorem but good numerical evidence may be considered sufficient). 

2. Without destroying well-posedness, modify your equations and choose your 
gauges, such that your problem actually becomes well-conditioned [J. 

3. Construct a solid numerical implementation, flexible enough to handle ex- 
periments as required by science and by finding solutions to the problems 
associated with point two. 

4. Discover (new) results in physics. 

5. Explain what you achieved (and how) to fellow numerical relativists and oth- 
ers, such as mathematical relativists, astrophysicists, cosmologists, or math- 
ematicians. 

Even without considering the last point (which the present article humbly tries 
to serve), numerical relativity is a challenging enterprise. 

The conformal approach complies with point one in the optimal sense: the 
equations are regular in the whole spacetime, including the asymptotic region, 
there are no ambiguities associated with ad-hoc cutoffs at finite distance, and the 
evolution equations are symmetric hyperbolic, which guarantees well-posedness 
of the initial value problem and allows to obtain well-posed initial- value-boundary 
problems. 

Point two however already poses a significant challenge: Well-defined is not 
well-conditioned, well-posed problems may still be hopelessly ill-conditioned for 
numerical simulation. A simple example is provided by any chaotic dynamical 
system (in the sense of ordinary differential equations). When it comes to solv- 
ing the Einstein equations, the gauge freedom of the theory results in having 
more equations (constraints and evolution equations) than variables, and more 

5 111 conditioned problems are those where a result depends very strongly on input, 
i.e. on initial data, see e.g. Sec. 1.6 and 6.1 of 152]. 
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variables than physical degrees of freedom. This redundancy can easily lead to 
spurious approximate solutions. Different ways to write the equations are only 
equivalent with regard to exact solutions, but approximations will tend to exhibit 
constraint violating or gauge modes that may grow very fast (e.g. exponentially). 
This is perfectly consistent with well-posedness but not acceptable numerically. 
Even without triggering instabilities, the choice of a bad gauge is likely to create 
features in the solution which are in practice impossible to resolve. The "good 
news" is that many of the problems encountered with the conformal field equa- 
tions have counterparts in traditional approaches to numerical relativity. The 
way toward solving these problems usually takes the form of gaining insight 
from simplifications and analytical studies, which then have to be tested in full 
numerical simulations. This requires a flexible code that is geared toward per- 
forming the necessary experiments, which leads to point three - another hard 
task for classical relativists, because it requires an engineering attitude many 
relativists are not familiar with. The gauge freedom of general relativity and 
absence of a natural background creates an additional twist when it comes to 
point four, which leads to numerous technical and conceptual subtleties. 

What is the roadmap for the future? In order to comply with points two 
and three of the above recipe, preliminary work is carried out toward a new 3D 
code that will be flexible enough to carry out a range of numerical experiments 
in order to come up with well-conditioned algorithms for the conformal field 
equations. One major issue in the improvement of algorithms is to implement a 
better boundary condition, that does not require a transition zone, allows the 
boundary to be closer to ^ + and minimize constraint violations generated at the 
boundary or outside J? . Here an essential problem is that has spherical cuts, 
and algorithms based on Cartesian grids are probably not optimal. Certainly, 
a lot of energy will have to be devoted to the question of finding appropriate 
gauge conditions. Particularly hard seems to be the question of how to choose 
the Ricci scalar R of the unphysical spacetime. Since R steers the conformal 
factor implicitly through nonlinear PDEs, it seems very hard to influence the 
conformal factor in any desired way. 

An important role in improving the analytical understanding and in set- 
ting up numerical experiments will be played by the utilization of simplifica- 
tions. Particularly important are spacetime symmetries and perturbative stud- 
ies. Minkowski and Kruskal spacetimes provide particularly important cases to 
be studied in this context. An alternative route to simplification, which has 
been very successful in numerical relativity, is perturbative analysis, e.g. with 
Minkowski or Schwarzschild backgrounds. In the context of compactification this 
has been carried out numerically with characteristic codes in [|53|]54|| (using ap- 
propriate variables in the Teukolsky equations, the perturbation equations are 
made regular at ^ + ), some of the problems that showed up there are likely to 
be relevant also for the conformal approach. 

The theory of general relativity is known as a never drying out source for 
subtle questions in physics and mathematics. Numerical relativity is hoped to 
help answer some important questions - but at the same time poses many new 
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ones. Without a thorough understanding of how to obtain approximate solutions, 
our insight into the theory seems incomplete. For isolated systems, the mastering 
of compactification techniques promises reliability and precision. The next years 
are expected to see some significant progress in this direction. 
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